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Section  1 


INTRODUCTION 

The  retrieval  of  atmospheric  temperature  and  moisture  profiles 
has  been  routinely  performed  using  data  from  orbiting  meteorological 
satellites  in  the  past  15  years.  However,  the  techniques  and 
procedures  used  in  the  retrieval  were  mostly  statistical  in  nature  and 
were  not  based  on  the  principles  of  radiative  transfer.  The  solution 
of  the  radiative  transfer  equation  for  the  Planck  function  distribution 
given  a  set  of  upwelling  radiances  at  the  top  of  the  atmosphere  is 
governed  by  a  Fredholm  equation  of  the  first  kind.  This  equation  is 
known  to  be  mathematically  ill-conditioned  because  of  the  instability 
introduced  by  the  inverse  operation. 

In  most  retrieval  approaches,  certain  a  priori  assumptions  or 
constraints  must  be  imposed  on  the  temperature  profiles  to  be 
retrieved.  For  example,  in  the  constrained  linear  inversion  method 
developed  by  Twomey  (1963)  ,  a  mean  temperature  profile  is  first 
assumed.  Subsequently,  the  retrieved  temperatures  are  obtained  using  a 
least-square  method  with  quadratic  constraints.  In  the  statistical 
method  (Liou,  1980),  the  difference  between  the  predicted  and 
climatological  mean  temperature  profiles  is  assumed  to  be  a  linear 
combination  of  the  deviation  of  measured  radiances  from  the  mean  for 
all  sounding  channels.  In  the  Backus -Gilbert  inversion  method 
(Conrath,  1972),  certain  assumptions  are  imposed  on  the  shape  of  the 
kernel  function.  In  other  numerical  methods  (Chahine ,  1970;  Smith, 
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1970) ,  a  simple  relationship  between  the  measured  radiances  and 
temperatures  is  assumed. 

In  an  attempt  to  break  away  from  the  conventional  techniques 
involving  profile  fitting  with  prechosen  sets  and  regression  methods 
based  on  statistics,  King  (1985)  proposed  a  novel  approach  to  solve  the 
temperature  inversion  problem.  Unlike  the  above-mentioned  methods,  the 
approach,  referred  to  as  the  differential  inversion  method  (DIM), 
directly  solves  for  temperatures  at  a  number  of  pressure  levels  that 
correspond  to  the  peaks  of  the  weighting  function.  The  DIM  is  an  exact 
method  for  the  inversion  of  the  temperature  profile,  in  which  the  local 
Planck  Intensity  can  be  exactly  expressed  by  a  linear  combination  of 
the  derivatives  of  upwelling  radiances  in  the  logarithmic  pressure 
coordinate.  The  method  is  free  from  the  need  for  a  priori  data  basing 
and  the  inversion  requires  no  constraints. 

In  this  report,  we  explore  the  generalization  and  practicality 
of  the  DIM  for  temperature  retrievals.  In  Section  2,  we  describe  the 
fundamentals  of  the  DIM.  In  Section  3,  applications  of  this  method  to 
synthetic  temperature  retrievals  using  HIRS  channels  are  illustrated. 
Finally,  conclusions  are  given  in  Section  4. 


2 


Section  2 


THE  PHYSICAL  FUNDAMENTALS  OF  THE 
DIFFERENTIAL  INVERSION  METHOD 
2 . I  Forward  Radiative  Transfer  Equation 

The  upwelling  radiance  at  the  top  of  the  atmosphere,  Ru,  for  a  given 
channel  may  be  derived  from  the  basic  radiative  transfer  equation  for  a 
plane-parallel  atmosphere  ur.^ier  local  thermodynamic  equilibrium.  In  the 
pressure  coordinate,  we  have 

R„  -  *„(TS)  r„(ps)  +  f°  Bu( p)  aTt/(p)  dp  ,  (2.1) 

JPs  dp 


where  B u  is  the  Planck  intensity  at  wavenumber  v,  T v  the  transmittance,  ps 
the  surface  pressure,  and  Ts  the  surface  temperature.  The  first  and  second 
terms  represent  surface  and  atmospheric  emission  contributions,  respectively. 

The  layer  below  the  surface  may  be  viewed  as  an  infinite  isothermal 
emitter  with  a  temperature  Ts .  Thus,  the  surface  emission  contribution  to 
the  radiance  R v,  as  represented  by  the  first  term  in  Eq.  (2.1),  may  be 
expressed  as 


5„(TS>  T„(ps) 


B„(p)  3Tt/(p)  dp 

dp 


(2.2) 


In  presenting  Eq .  (2.2),  it  is  noted  that  Tu(^)  -  0.  Also,  Bu{ p)  ~  5^(TS) 
for  ps  <  p  <  <».  In  the  band  center,  T„(ps)  -  0,  so  that  surface  contribu¬ 
tion  to  the  randiance  R„  may  be  neglected.  However,  in  the  wing  of  a  band, 


7\,  ( ps)  >  0,  the  surface  emission  could  be  an  important  source  for  the 
upvelling  radiance.  In  either  case,  we  may  substitute  Eq .  (2.2)  into  Eq . 
(2.1)  to  obtain 


Ru 


Bu(  p) 


ar„(P) 

dp 


dp 


(2.3) 


The  weighting  function,  which  signifies  the  weight  of  the  Planck 
intensity  contribution  to  the  upwelling  radiance,  may  be  defined  in  the 
logarithm  of  pressure  in  the  form 


V„(p) 


ar„(p) 

dinp 


(2.4) 


For  each  wavenumber  u ,  the  transmittance  Tv{ p)  decreases  from  1  to  0  as  the 
parameter  inp  increases  from  -<*>  to  «.  Hence,  the  weighting  function  is 
always  a  positive  quantity,  and  there  must  be  a  peak  value  located  at  p  -  p , 
where  the  rate  of  change  of  Ti/( p)  reaches  a  maximum.  The  peak  pressure 
level  p  varies  with  wavenumber  v.  If  v  is  in  the  band  center,  p  will  be 
located  at  a  relatively  low  pressure  level,  whereas  if  v  is  in  the  wing 
region,  p  will  be  near  the  surface. 

In  practice,  a  satellite  sounding  instrument  spans  a  spectral 
interval.  The  sensitivity  of  the  instrument  varies  with  the  wavenumber, 
according  to  a  response  function  <j> .  Assuming  that  there  are  M  channels, 
the  radiance  measured  by  the  ith  channel,  with  response  function 
may  be  expressed  as  follows: 


R 


i-r 

J  n 


B  i(p)  f/i(p)  dp/p 


i  -  1,2, 


,M, 


(2.5) 
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where  Che  weighting  function  W^(p)  is  now  the  average  of  W„(p) ,  defined  in 
Eq .  (2.4),  weighted  by  the  response  function  ,  viz. 

^i(p)  “  f  4LCv)  dv  ,  (2.6) 

J  Ai/j_ 

with  Aisi  the  spectral  interval  for  the  ith  channel.  In  writing  Eq .  (2.5), 
we  omit  the  wavenumber  dependence  of  the  Planck  intensity  since,  over  a 
small  spectral  interval,  the  Planck  intensity  does  not  vary  significantly 
with  frequency.  Equation  (2.5)  is  a  Fredholm  equation  of  the  first  kind. 

The  conventional  approach  for  solving  this  kind  of  equation  is  to  assume  a 
functional  form  for  Bj_(p),  and  then  determine  the  unknown  coefficients  in 
the  function,  based  on  the  given  distributions  of  R^  and  In  practice, 

since  only  a  finite  number  of  R^  values  are  available,  the  solution  of  B^( p) 
from  the  forward  radiative  transfer  equation  is  mathematically  ill- 
conditioned.  However,  if  the  integration  can  be  removed  by  a  proper  inverse 
transformation,  then  a  priori  assumptions  or  constraints  commonly  added  to 
the  retrieved  temperature  profile  will  no  longer  be  required. 

2 • 2  Differential  Inversion  Method 

We  shall  approach  the  inverse  problem  by  a  transformation  of  var¬ 
iables.  We  first  define  the  scaled  pressure  p  -  p/p^,  Since  the 
weighting  function  W^(p)  peaks  at  p^(p  -  1) ,  we  may  write  W^(p)  -  max 
[Wi(p)].  Furthermore,  the  following  variables  are  introduced: 

n  -  -  inp  , 
n  -  -  inp  , 

v  -  7T  -  it  —  in(p/p)  (2.7) 
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In  terms  of  these  new  variables,  Eq .  (2.5)  may  be  rewritten  in  the  form 


( 7r )  -  J°°  B(ir)  W(jt-it)  d?r 


(2.8) 


where 


R(*r)  -  «i(p) 
B(ir)  -  fli(p) 
W(-w)  -  W(p) 


'.'e  then  perform  a  bilateral  Laplace  transformation  on  both  sides  of  Eq 
(2.8),  viz . , 


(2.9) 


e  S?r  R(jt)  dff 


{L  B  ( 7r )  W(jr-rr)  d*  ^  dn 


(2.10) 


By  virtue  of  the  convolution  theory  for  Laplace  transform  (Widder,  1971),  we 
may  deconvolute  the  right-hand  side  of  Eq.  (2.10)  to  obtain 


r(s)  -  b(s)  w(-s) 


r(s)  -  e  Sn  R ( 7r )  dw 


(2.11) 


(2.12) 


e  Sn  B  ( 7T )  d7T 


(2.13) 


w(— s) 


(2.14) 


*CO 

eSn  W(w)  dn 

— ao 

In  deriving  Eq.  (2.11),  there  are  a  number  of  intricate  sign  changes  so  that 
the  argument  for  w  is  -s  rather  than  s,  and  all  the  integrations 
are  over  the  w-space.  The  detailed  derivations  are  in  Appendix  A.  The 
quantities  r(s),  b(s),  and  w(-s)  are  the  bilateral  Laplace  transforms  of  R, 
B,  W,  respectively.  Since  the  integration  of  B  and  W  in  Eq.  (2.10)  has  been 
deconvoluted,  we  may  express  b(s),  which  is  related  to  the  Planck  intensity, 
in  the  form 


b(s) 


r(s)_ 
w(— s) 


(2.15) 


so  that 


B(ir) 


r(s) 

w(-s) 


(2.16) 


where  L  ^  is  the  inverse  bilateral  Laplace  transform.  It  is  clear  that  once 
the  functional  forms  for  r(s)  and  w(-s) ,  as  well  as  the  inverse  formula  for 
r(s)/w(-s)  are  known,  the  profile  of  B(7r)  may  be  derived.  One  approach 
would  be  to  expand  the  quantity  l/w(-s)  into  a  MacLaurin  series  (Pearson, 
1974),  which  is  a  form  of  Taylor’s  expansion  with  the  reference  point  fixed 
at  s  -  0,  in  the  form 


1 

w(-s) 


CO 


z 

k-0 


*k  sk 


(2.17) 


where  the  coefficient  is  related  to  the  kth  derivative  of  the  function 
l/w(-s)  evaluated  at  s  -  0,  as  follows: 
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i 


1  1 


w(-s) 


(2.18) 


Equation  (2.17)  may  be  substituted  into  Eq .  (2.16)  to  obtain  an  expression 
for  b(s)  in  terms  of  an  infinite  power  series  in  the  form 


b  ( s )  -  £  Ak  sk  r(s) 

k-0 


(2.19) 


Thus,  from  Eq .  (2.16),  we  have 


—  -1 
B(w)  “  X  Ak  L  [sk  r (s) 


The  Laplace  transform  of  the  kth  derivative  of  a  function  R(x)  can  be 
expressed  by  (Abramowitz  and  Stegun,  1968) 


"  dkR(x)  1  k-1  .  , 

L  ~TTk  -  sk  r(s)  -  l  sn  R(k~:  n)(x  -  — ) 

L  J  n-0 


(2.20) 


(2.21) 


when  7r  -+  -®)  i.e.,  p  -♦  «,  the  upwelling  radiance  R (. 7r )  does  not  exist  in 
physical  space.  Thus,  from  Eq .  (2.21),  we  have 


L  ^s1"  r(s)  ]  -  R(k)(7r) 


(2.22) 


The  Planck  intensity  at  a  given  level  n  is  now  expressible  in  terms  of  the 
linear  sum  of  radiance  derivatives  at  that  level  in  the  form 


B(»)  -  l  Ak  R(k)  (*) 


(2.23) 
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The  solution  of  B(jt)  will  be  mathematically  exact  for  infinite  summations. 

It  requires  no  constraint  and  is  self-limited  to  the  highest  order  of 
recoverable  derivatives.  In  practice,  since  there  are  only  a  limited  number 
of  radiance  measurements  available  for  distinct  values  of  n,  two  error 
components  may  be  entered  into  the  retrieval.  One  is  the  instrument  noise, 

while  the  other  is  the  so-called  "null-space"  error,  which  is  generated  by 

00  ~  - 

an  inadequate  representation  of  R  (tt)  in  the  data-void  m-space.  Both  of 
these  errors  may  be  amplified  or  suppressed  by  the  coefficients  A^,  depending 
on  their  magnitude. 


2.3.  Weighting  Function  and  Inverse  Coefficient 

The  determination  of  B(jt)  from  Eq .  (2.23),  requires  knowledge  of  both 

(k)  - 

the  inverse  coefficient  A^  and  the  derivatives  of  radiance  R  (n) .  In  this 
section,  we  shall  focus  on  the  discussion  of  the  derivation  of  expressions 
for  A^.  We  shall  rewrite  Eq .  (2.14)  in  the  form 


w(-s) 


*00 

p  S  W(p)  dp/p 

0 


where  the  spectral  weighting  function  W(p)  is  defined  by 


(2.24) 


W(p)  - 

Note  that  in  Eq .  (2.24),  the  integration  is  over  the  p-space  rather 
than  the  w-space.  Since  w(-s)  is  a  function  of  s  only,  use  of  either  p 
or  n  is  mathematically  correct.  In  physical  space,  however,  it  is  more 
appropriate  to  express  W  in  terms  of  p.  In  principle,  if  the  weighting 


ar„(p) 


Sfnp 


4>i(v)  di/ 


(2.25) 
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functions  !7(p)  are  known,  e.g.,  through  accurate  line-by-line  spectral 
integration,  it  is  possible  to  evaluate  l/w(-s)  and  its  derivatives  numeri¬ 
cally.  Subsequently,  A^  can  be  computed  from  Eq .  (2.18).  In  this  manner, 
the  potential  error  in  A^  would  be  minimal.  However,  for  the  purpose  of 
analysis,  it  is  desirable  to  develop  an  analytic  form  that  can  approximate 
the  weighting  functions  associated  with  sounding  channels  without  incurring 
significant  errors.  A  generalized  weighting  function  was  proposed  by  King 
y 1985)  in  the  form 

W(p)  -  mm_1  r"1 <m)  p  exp  [-m  p  1/m ]  ,  (2.26) 

'.-.’here  T  is  the  Gamma  function  and  m  is  an  index  controlling  the  sharpness  of 

the  weighting  function.  When  m  -  1  and  0.5,  the  weighting  functions  follow 

Goody-statistical  and  Elsasser-regular  band  models,  respectively.  The 
m-1  —1 

factor  m‘  T  (m)  in  Eq.  (2.26)  normalizes  the  weighting  function  to  unity 

and  insures  that  W  -  W  -  W(p  -  1) . 

max  v 

Figure  1  depicts  the  shapes  of  the  generalized  weighting  functions  for  m 
“0.5,  1,  and  2.  Large  m  values  correspond  to  broader  weighting  functions, 
but  smaller  peak  values.  The  peak  values  for  m  -  0.5,  1,  and  2  are  0.484, 
0.368,  and  0.271,  respectively.  As  shown  in  Fig.  1,  the  prime  contribution 
of  the  weighting  function  comes  from  0.1  <  (p/p)  <  10,  a  range  of  two 
orders  of  magnitude. 

Using  the  generalized  weighting  function,  we  can  show  from  Eq .  (2.24) 

tha  t 

w(— s)  -  r(m(l-s)  )  /  [r(mfm_niS]  .  (2.27) 
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I' i t o  derivation  of  this  form  is  given  in  Appendix  B.  Evaluation  of  the 
inversion  coefficient  A^,  based  on  Eq .  (2.27),  requires  derivatives  of  the 

Camma  functions,  which  are  related  to  digamma  and  polygamma  functions.  The 
expressions  for  for  k  -  1  to  5  as  functions  of  the  index  m  have  been 
derived  and  are  given  below: 


12 


Tables  for  Che  values  of  0^)  (m)  are  available  for  1  <  m  <  2  and  k  <  3 
(Abramowitz  and  Stegun,  1968).  However  for  m<l,  m  >  2,  or  k  >  3,  these 
polygamma  functions  must  be  evaluated  using  recurrence  relationships  and 
series  expansions.  We  find 


,(k).  .  ,(k).  , .  .  . . k+1  ,  -k-1 

0  (m)  -  0  (m+ !)+(-!)  k!m 


m  <  1 , 


0(k)(m)  -  0(k)(m~l)  +  (-l)k  k!  (m-1)  k  1  ,  m  >  1, 


(2.34) 

(2.35) 


,  (k)  ,  , .  .  , . k+1 

0  '(m+1)  -  (-1) 


k!  f (k+1) 


■°7pj'  f(k+2)  m  +  r (k+3)  m2  -  ... 


,  k  >  3 


(2.36) 


where  f(k)  is  the  Riemann  Zeta  function,  which  is  defined  by 


r(k)  -  l  n  k  .  (2.37) 

n-1 

Figure  2  shows  the  inversion  coefficients  A^  (k  -  1-5)  as  functions  of 

the  index  m.  For  any  given  m,  may  be  obtained  from  this  figure.  For  0.1 

<  m  <  1,  A^  are  within  -1  and  0.2.  Based  on  Eq.  (2.23),  any  error  incurred 
(k-V 

in  R  (w)  will  be  substantially  reduced.  This  implies  that  for  0.1  <  m  < 

^  ’  —  He)  — 

B  ( 7r )  will  be  relatively  insensitive  to  errors  in  Rv  '(tt)  .  For  m  >  1,  |A^|, 

| A 3 | ,  and  | A 5 |  are  less  than  1.  However,  | A2 |  and  | A4 1  are  larger  than  1. 
Errors  in  R^(7r)  and  R^(jt)  will  then  be  greatly  amplified  through  A2  and 
A/(  and  propagate  into  the  computation  of  B(w) . 


13 


Section  3 


APPLICATIONS  OF  THE  DIM  TO  TEMPERATURE 
RETRIEVALS  USING  THE  HIRS  CHANNELS 

In  this  section,  we  shall  test  the  DIM  by  performing  synthetic 
analyses  using  HIRS  channels  in  the  15  pm  CO2  band.  It  is  our  objective 
to  physically  understand  the  generalization  and  determine  the  practi¬ 
cality  of  this  method  for  temperature  retrievals. 

3 . 1  HIRS  Channel  Characteristics 

In  order  to  investigate  the  potential  applicability  of  the  DIM,  HIRS 
channels  (Smith  et  al. ,  1975)  are  used  to  perform  the  temperature 
retrieval.  HIRS  consists  of  seven  channels  in  the  15  pm  CO2  band,  and 
five  channels  in  the  4.3  (im  CO2  band.  We  shall  focus  our  investigation 
on  the  seven  channels  in  the  15  ^m  CO2  band.  It  Is  noted  that  the 
peaks  of  the  weighting  functions  of  these  channels  are  more  evenly 
distributed  than  those  in  the  4.3  jun  band.  The  characteristics  of 
these  channels  are  listed  in  Table  1. 


Table  1.  HIRS  channel  characteristics. 


Channe 1 

v  (cm-^) 

v2 

A  v 

Principal 

Absorbers 

Level  of 
^max 

1 

668 

666 

670 

4 

C02 

30 

2 

679 

674 

684 

10 

C02 

60 

3 

690 

685 

697 

12 

C02 

100 

4 

702 

696 

712 

16 

C02 

250 

5 

716 

708 

724 

16 

CO  2 

500 

6 

732 

724 

740 

16 

co2/h2o 

750 

7 

748 

740 

756 

16 

co2/h2o 

900 

Computation  of  Weighting  Functions  by  the  K-Correlated  Method 


In  order  to  compute  the  transmittances  and  weighting  functions,  we 
used  the  absorption  coefficient  datasets  computed  by  Chou  and  Kouvaris 
(1986)  based  on  line-by-line  data  compiled  by  Rothman  et  al .  (1983). 

The  absorption  coefficients  due  to  CO2  absorption  are  available  for 
every  0.002  cm- ^  between  540  and  800  cm-^  at  19  pressure  levels  (0.25- 
1000  mb  for  every  Alogp  -  0.2)  and  three  reference  temperatures  (210, 
250,  and  290  K) .  For  each  wavenumber  and  pressure  level,  the  computed 
absorption  coefficients  were  fitted  to  an  empirical  form  of  the  tempera 
ture  adjustment  as  follows: 

k„  -  exp  [a  +  bAT  +  cAT^]  ,  (3.1) 

where  AT  -  T-250.  The  coefficients  a,  b,  and  c,  obtained  by  curve 
fitting,  are  pressure  and  wavenumber  dependent.  Thus,  the  absorption 
coefficient  for  any  wavenumber,  pressure,  or  temperature  can  be 
obtained  by  interpolation  between  reference  values. 

Using  the  preceding  absorption  coefficients,  the  spectral  transmit¬ 
tance  and  weighting  function  for  a  given  pressure  level  p^  and  for  each 
HIRS  channel  may  be  computed  by  the  following  procedures.  First,  a 
table  of  k-distribution  values  for  each  of  the  19  pressure  levels  and 
for  each  channel  is  compiled  based  on  a  reference  atmospheric  profile. 
The  k-distribution  is  basically  the  probability  of  a  certain  k  value  in 
a  small  wavenumber  interval.  The  determination  of  spectral  transmit¬ 
tance  over  this  interval  does  not  depend  on  the  position  of  each  k 
value,  but  on  its  frequency  of  occurrence  (Arking  and  Grossman,  1972). 


Let  f(k)  be  the  probability  density  function,  we  may  define  a  cumulative 
probability  density  function  in  the  form 


g(k) 


■k 

f(k')  dk' 

0 


(3.2) 


It  is  noted  that  g(0)  -  0  and  g(»)  -  1.  In  practice,  for  each  channel 
and  pressure  level,  we  divide  the  logk-space  into  a  large  number  of 
intervals  of  equal  width.  Using  the  k  values  calculated  for  every 
0.002  cm-l,  we  calculate,  in  a  discrete  manner,  the  k-distribution 
function  for  the  ith  Alog  k  interval  in  the  form 


f(kL)  Alog  ki  -  n(ki)  ,  (3.3) 

N 

where  n(k*)  is  the  number  of  wavenumber  intervals  (Ai/  -  0.002  cm-'-) 
whose  logk  value  falls  between  logk^  and  logk^  +  Alogk^,  and  N  is  the 
total  number  of  wavenumber  intervals  for  a  given  channel.  The  cumula¬ 
tive  g-functioi.  ip  to  the  ith  interval  can  be  obtained  from 

g(ki)  -  l  __J_  ,  i  >  J  •  (3. A) 

j-1  N 

The  next  step  for  calculating  the  transmittance  or  weighting 
function  is  based  on  the  correlated  k-distribution  method  (Lacis  et 
al.,  1979).  For  a  homogeneous  path,  the  spectrally  averaged  transmit¬ 
tance  may  be  expressed  by 

r_< u)  -  f  e-k*u  f-  ■  <35) 

V  J  .  Lv 

J  Al/ 
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As  mentioned  above,  over  a  small  spectral  interval  where  the  Planck 
intensity  may  be  considered  constant,  the  position  of  each  k  value  has 
no  bearing  on  the  calculation  of  the  spectral  transmittance.  Thus,  the 
spectral  integral  may  be  transformed  into  the  integral  over  the  k- 
domain  as  follows: 


rk. 


r_(u)  - 


max 


^min 


e  ku  f(k)  dk 


(3.6) 


where  f(k)  -  (dk/di/)  ^/t\v.  Based  on  the  definition  of  the  g-function 
denoted  in  Eq.  (3.2),  Eq .  (3.6)  may  be  rewritten  in  the  form 


T_(u) 

v 


*1 

0 


e~k(g)u 


dg 


(3.7) 


If  the  relationship  between  k  and  g  is  known,  the  spectral  transmittance 
T_(u)  can  be  evaluated. 

However,  for  an  inhomogeneous  atmosphere,  two  assumptions  must  be 
made  in  order  to  obtain  a  form  similar  to  Eq.  (3.7).  The  spectral 
transmittance  for  an  inhomogeneous  atmosphere  may  be  expressed  by 


r_(u) 


-  exp  -  k^  (p ,  T) 

J  Ai/  L  JO 


du' 


di/ 

Ai/ 


(3.8) 


•  Let  pr  and  Tr  denote  the  reference  pressure  and  temperature,  respective¬ 
ly.  Then  for  any  pair  of  (I'i.i'j),  if  k„(pr , Tr)  Ji/^  -  k„(pr ,Tr)  |i/j  ,  it  is 
assumed  that  k^Cp.T)^^  -  k„(p  ,T)  |«/j  .  This  assumption  implies  that  when 

•  the  pressure  and  temperature  change  from  the  reference  values  (pr,Tr) 
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to  arbitrary  values  (p,T),  the  ratio  k(kr,P,T)  -  k^p.T)  /  k^Cpj-.Tj.)  is 
a  function  of  k^ (pr , Tr) («kr) ,  p  and  T.  Since  kr  are  known  values,  fc  is 
independent  of  the  wavenumber.  Thus,  the  pathlength  integral  in  the 
exponential  argument  in  Eq.  (3.8)  may  be  scaled  in  the  form 

u 

Mp.T)  du'  -  kr  {i  <kr,p,T)  ,  (3.9) 

.  0 

where  the  scaled  path  length  u  is  defined  by 

*u 

u(kr,p,T)  -  £  (kr,p,T)  du'  .  (3.10) 

.  0 

In  this  way,  the  wavenumber  dependence  in  kr  is  separated  from  the 
pressure  and  temperature  dependence  in  the  pathlength  integration.  We 
may  then  transform  the  spectral  integral  in  Eq.  (3.8)  into  a  k-space 
integral  to  obtain 

T_(u)  -  exp  [~k(gr)  u]  dgr  ,  (3.11) 

where  gr(k)  is  the  accumulated  probability  density  function  for 
reference  pressure  and  temperature.  To  perform  numerical  calculations 
using  Eq.  (3.11)  further  assumption  must  be  made.  For  any  pair  of  (i/^, 
uj),  if  ky  (Pr.TjOli/i  >  k„(pr  ,Tr)  |  i/j  ,  then  k^  ( p ,  T )  |  ^  >  k^(p  ,  T)  |  i/j  . 

This  assumption,  together  with  the  previous  one,  allows  us  to  express 
the  following  relationship: 


fr(k' )  dk'  -  J*  f(k' )  dk' 


or ,  from  Eq .  (3.2), 


gr(k)  -  g(k) 


(3.12) 


(3.13) 


Based  on  che  preceding  assumptions,  Eq.  (3.11)  may  be  rewritten  in  the 
form  similar  to  Eq .  (3.7)  for  homogeneous  atmospheres: 


r_(p)  -  exp  [-k(g)  u]  dg(k) 

*  Jo 


(3.14) 


In  practice,  we  use  the  g-function  table  for  various  pressure  levels 
compiled  in  the  first  step,  along  with  Eq.  (3.13),  to  determine  the 
correspondence  between  k  and  kr.  Then,  for  each  g  value,  we  determine 
the  kr  value  and  a  set  of  £(kr,p,T)  to  compute  the  pathlength  inte¬ 
gration.  Finally,  the  spectral  transmittance  is  obtained  by  numerical 
integration  according  Eq.  (3.14).  This  is  referred  to  as  the  correlated 
k-distribution  method.  Lacis  et  al.  (1979)  showed  that  this  method  is 
remarkably  accurate  for  the  computation  of  transmittances  and  heating 
rates  for  the  9 . 6  ^ra  ozone  band. 

In  the  present  synthetic  calculation,  the  transmittance  for  each  HIRS 
channel  is  computed  according  to  Eq.  (3.14).  The  transmittance 
calculations  using  the  correlated  k-distribution  are  extremely  efficient 
and  accurate.  The  results  are  almost  identical  to  those  from  line-by- 
line  calculations  (Liou  et  al.,  1988).  Figure  3  shows  the  weighting 
functions  for  seven  HIRS  channels,  based  on  the  U.S.  standard  atmos¬ 
pheric  profile.  The  peak  levels  of  these  weighting  functions  are 


exactly  the  same  as  those  listed  In  Table  1.  It  is  noted  that  channel 
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« 


1  peaks  at  the  highest  level  and  has  the  broadest  weighting  function, 
whereas  channel  7  peaks  near  the  surface  and  has  the  most  narrow 


I 

weighting  function. 

3 . 3  Computation  of  Inverse  Coefficients  Using  the  Generalized 
Weighting  Function 

It  is  possible  to  compute  A^  directly  from  Eqs .  (2.14)  and  (2.18) 
using  the  weighting  functions  obtained  from  the  correlated  k-distribu- 
tion  method,  but  the  computation  could  be  very  tedious.  With  the  use 
of  the  generalized  weighting  function  defined  in  Eq .  (2.26),  however, 

A^  can  be  evaluated  according  to  Eqs.  ( 2 . 28 ) - ( 2 . 32 )  in  a  straightforward 
manner.  We  have  developed  a  numerical  method  to  fit  HIRS  weighting 
functions  to  the  generalized  form  to  obtain  the  index  m  associated  with 
each  channel.  It  is  similar  to  the  least-square  method.  We  define  the 
sum  of  the  weighted  square  error  at  a  given  pressure  level  in  the 
form 


E 


^m<Pi) 


W(Pi)]2  «i 


(3.15) 


where  L  is  the  total  number  of  pressure  levels,  Pi  -  p^/p,  and  W 
are  the  generalized  and  computed  weighting  functions,  respectively. 

The  weight  factor  «i  is  prescribed  as  follows: 

f  exp  (-p£)  ,  Pi  >  1 

n-\  -i  .  (3. if) 

l  exp  (-pi  )  ,  Pi  <  1 


In  this  way,  more  weight  is  placed  on  the  error  near  the  weighting 
function  peak  (pi  -  1)  than  in  the  wing  region.  This  is  appropriate 


since,  in  computing  the  upwelling  radiance,  the  greatest  contribution 
comes  from  the  weighting  function  peak  level.  We  then  proceed  to 
minimize  E,  such  that 


3E 

3m 


-  0  , 


(3.17) 


which  leads  to 

L 

Z  2  [ ^ rn ( P i ^  -  ^(Pi)l  £i  ^m(Pi) 

i-1 

(1  -  1/m  -  (pL)1/m  [1  -  in  (pi)/m]  +  inm  -  0(m) )  -  0,  (3.18) 

where  t/>(  m)  is  the  digamma  function.  The  derivation  of  Eq .  (3.18)  is 
given  in  Appendix  C. 

Equation  (3.18)  is  a  complicated,  nonlinear  algebraic  equation,  and 
can  only  be  solved  by  a  numerical  root-finding  method.  Although 
Newton's  iteration  method  is  the  most  efficient,  it  requires  specifica¬ 
tion  of  the  functional  form  of  the  second-order  derivative  of  E,  which 
in  the  present  case  is  much  more  complicated  than  the  f.  jt-order 
derivative.  For  any  channel,  m  will  be  within  0.1  and  4  as  shown  in 
Table  2.  Thus,  we  may  use  the  bisection  method  described  by  Hamming 
(1973)  to  solve  for  m.  The  details  of  the  bisection  method  used  in  the 
fitting  are  given  in  Appendix  D,  Table  2  lists  the  values  of  the 
parameter  m  for  each  HIRS  channel,  RMS  error,  peak  value  of  f/(pj_),  and 
error  near  the  peak  of  the  weighting  function.  The  RMS  error  is 


defined  by 


Table  2. 


Values  of  parameter  m,  RMS  error,  and  error  near  the  peak 
weighting  function  for  the  HIRS  channels. 


Channel 

m 

RMS  error 

^max* 

error  near  peak 

1 

2 . 8370 

0.0141 

0.234 

-0.0036 

2 

0.6410 

0.0137 

0.433 

0.0073 

3 

0.6668 

0.0107 

0.436 

0.0017 

4 

0.4570 

0.0237 

0.487 

-0.0129 

5 

0.4273 

0.0137 

0.517 

-0.0048 

6 

0.2305 

0.0159 

0.614 

0.0076 

7 

0.3160 

0.0118 

0.566 

-0.0004 

*{Vmax  denotes 

the  maximum 

value  of  the 

weighting 

function . 

,  L  1/2 

e“  7  l  [Wm(pi)  -  W(Pi)]2  -  (3.19) 

1  L  i-1 

The  RMS  errors  for  all  channels  are  within  about  0.01-0.02.  Such  errors 
are  sufficiently  small  to  have  a  significant  effect  on  the  upwelling 
radiance  calculation.  The  effect  of  the  weight  factor  is  reflected 
by  the  fact  that  the  error  near  the  peak  for  each  channel  is  smaller  in 
magnitude  than  the  corresponding  RMS  error.  For  channel  1,  a  large  m 
is  found,  which  corresponds  to  a  broad  weighting  function.  For  other 
channels,  m  ranges  from  0.2  to  0.7.  It  is  apparent  that  neither  the 
random  model  nor  the  regular  model  can  fit  the  HIRS  channel  transmit¬ 
tance  well.  Based  on  the  m  values,  the  inversion  coefficient  A^  may  be 
evaluated  for  each  channel  (see  Fig.  2).  These  A^  values  are  tabulated 
in  Table  3.  Note  that  all  |AjJ  are  less  than  1,  except  | A2 |  and  | A^ | 
for  channel  1.  Thus,  it  is  anticipated  that  the  retrieved  temperature 
at  30  mb,  which  is  the  peak  level  of  weighting  function  for  channel  1, 
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Table  3. 


Inversion  coefficient  values,  for  the  seven  HIRS 

channels . 


I 


Channel 

A1 

a2 

a3 

a4 

a5 

1 

-0.52 

-1.60 

0.21 

1.44 

0.270 

2 

-0.60 

-0.48 

0.01 

0.06 

0.017 

3 

-0.61 

-0.50 

0.01 

0.07 

0.018 

4 

-0.63 

-0.39 

0.00 

0.04 

0.009 

5 

-0.63 

-0.37 

0.00 

0.04 

0.007 

6 

-0.70 

-0.27 

-0.01 

0.03 

0.001 

7 

-0.67 

-0.32 

-0.01 

0.03 

0.003 

will  be  very  sensitive  to  the  measurement  noise  and  errors  in  the 
radiance:  fitting. 

3 . 4  Synthetic  Computation  of  Channel  Radiances 

To  test  the  reliability  and  applicability  of  the  DIM,  real-time 
satellite  radiance  data  must  be  used.  At  the  same  time,  the  retrieved 
temperature  profiles  must  be  verified  with  the  ground- truth.  Before 
the  application  of  real  data,  we  shall  use  a  set  of  radiance  values 
obtained  from  synthetic  computations  to  explore  the  feasibility  of  the 
DIM  for  temperature  retrievals.  Let  R ^  be  the  radiance  associated  with 
the  ith  sounding  channel  whose  value  may  be  computed  from  Eq.  (2.5). 

The  U.S.  Standard  Atmospheric  Temperature  is  used  as  the  reference 
profile.  Table  4  lists  the  computed  upwelling  radiances  for  the  seven 
These  radiances  increase  from  the  band  center  to  the  wing. 


Table  4. 


Radiance  values  for  HIRS  channels. 


Channel 


Radiance 

R i  (erg/cm2-sec-sr) 


1 

224 

2 

437 

3 

513 

4 

818 

5 

999 

6 

1^30 

7 

1420 

3 . 5  Polynomial  Fitting  of  Radiances 

Tt  is  necessary  to  develop  a  curve- fitting  method  to  fit  the  seven 
upwelling  radiances,  since  higher-order  derivatives  are  required  to 
calculate  the  Planck  intensity.  Several  numerical  methods  are  available 
for  the  present  purpose,  and  can  be  divided  into  two  categories.  The 
first  is  the  spline  method,  which  fits  selected  radiance  values  in  a 
piece -wise  manner.  The  curve  segments  between  successive  radiance 
values  are  smooth  without  oscillations  incurred  by  higher-order 
polynomial  fittings.  If  the  cubic  spline  method  is  used,  an  assumption 
about  the  derivatives  at  the  endpoint  must  be  made.  We  may  obtain  at 
most  third-order  derivatives,  since  all  higher-order  derivatives 
vanish.  The  second  category  is  associated  with  the  polynomial  fitting. 
This  method  fits  all  radiance  values  to  a  single  polynomial. 

In  the  present  study,  we  use  the  Newtonian  interpolation  scheme 
(Hamming,  1973)  in  view  of  its  simplicity  in  the  determination  of 
higher-order  derivatives.  The  Newtonian  interpolation  scheme  uses  a 
polynomial  in  the  form 


(3.20) 


*  (k)k__ 

R  -  R\  +  L  II  (n  -  iri) 

k-1  j-1 


where  R  ]_  is  the  radiance  for  channel  1,  K  the  order  of  the  polynomial, 

(k)  - 

the  pseudo-derivative  of  order  k,  n  -  -  inp,  and  ttj  -  -  inpj  .  The 

(k) 

pseudo-derivatives  R are  defined  by 


(1)  r2  ~  Rl 

R  i  "  - 

*2  -  *1 


(1)  Ri+1  -  R  i 

Ri  ~  - 

*i+l  -  Iti 


1  <  i  <  K-1 


(2) 

*1 


(1)  (1) 

R2  -  R i 


*3  -  *1 


(1)  (1) 

(2)  *i+l  ~  *i 

Rl  -  _  ,  1  <  i  <  K-2 

jri+l  -  *£ 


(K-2)  (K-2) 

(K-1)  r2  ~  Rl 

Rl  -  _  ,  (3.21) 

*K  ~  *1 

where  Ri  is  the  radiance  for  the  ith  channel.  On  the  basis  of  the 
preceding  definition,  if  K  -  1,  the  polynomial  will  give  the  exact 
radiance  values  for  the  first  two  channels.  If  K  -  6,  the  exact 
radiance  values  for  all  seven  channels  will  be  produced. 

Figure  4  shows  the  curves  of  the  fitted  first-,  third-,  and  fifth- 
order  polynomials  based  on  the  seven  radiance  values  listed  in  Table  4. 
It  is  evident  that  the  fifth-order  polynomial  fits  the  radiance  values 
best.  However,  fluctuations  occur  between  channels  1  and  3.  This 
suggests  that  a  smooth  curve  -  fitting  method  could  be  advantageous  in 
the  inversion  exercises  because  higher-order  derivatives  of  the  curve 
are  needed  in  order  to  compute  the  Planck  intensity.  Based  on  Eq. 
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FIGURE  4.  The  computed  upwelling  radiances  (circles)  for  seven  HIRS 
channels  in  the  logarithmic  pressure  coordinate,  and  curve 
fittings  to  the  computed  values  using  first-,  third-,  and 
fifth-order  polynomials.  The  U.S.  standard  temperature  is 
used  in  the  computation. 


(3.20),  analytic  expressions  for  the  higher-order  derivatives  may  be 
derived  and  are  given  in  Appendix  E. 

3 . 6  Inversion  Exercise  Using  the  Generalized  Weighting  Function 

Since  the  inversion  coefficients  A^  are  known,  as  are  the  deriva¬ 
tives,  the  Planck  intensity  B(nj_)  at  the  peak  of  the  weighting  function 
for  each  channel  can  be  computed  from  Eq.  (2.23).  From  the  Planck 
intensity,  the  temperature  can  be  evaluated  from 

T(jri)  -  he  vi  /  {K  in  [1  +  2hcz  /  B  (nt) )  )  ,  (3.22) 

where  the  quantities  h,  c,  17^,  and  K  represent  the  Planck  constant, 
speed  of  light  in  vacuum,  central  wavenumber  of  the  ith  channel,  and 
Boltzmann  constant,  respectively. 

Figure  5  shows  retrieval  results  for  the  first-,  third-,  and  fifth- 
order  polynomial  fittings  to  synthetic  HIRS  radiances  using  the  U.S. 
standard  atmospheric  profile.  It  is  quite  encouraging  to  find  that  the 
retrieval  results  converge  to  the  true  temperature  in  the  lower 
troposphere  as  higher-order  radiance  derivatives  are  incorporated  in 
the  calculation.  For  the  fifth-order  approximation,  errors  in  the 
retrieved  temperatures  in  the  lower  troposphere,  corresponding  to 
channels  4  to  7 ,  are  within  about  2  K.  Retrieved  temperatures  associ¬ 
ated  with  channels  1-3  show  larger  deviations  from  the  true  temperature 
values.  The  large  deviations  in  the  upper  atmosphere  are  primarily 
caused  by  the  "null-space"  error  of  the  higher-order  derivatives, 
associated  with  the  less  satisfactory  polynomial  fitting  between 
radiances  for  channels  1  and  3,  as  noted  in  Fig.  4. 
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To  further  test  the  practical  applicability  of  the  DIM,  we  performed 
a  similar  retrieval  exercise  using  the  tropical  atmospheric  profile 
(McClatchey  et  al . ,  1971).  As  shown  in  Fig.  6,  the  tropical  temperature 
profile  differs  significantly  from  the  standard  atmospheric  profile. 

Most  notable  is  the  strong  inversion  of  the  tropopause  near  the  100  mb 
level.  Again,  the  retrieved  solution  converges  to  the  true  profile  in 
the  lower  troposphere  when  the  higher-order  derivatives  are  added  to 
the  series  expansion.  The  retrieved  temperatures  are  within  about  1  K 
for  channels  4-7.  However,  errors  in  the  temperature  at  the  top  three 
levels  are  even  larger  than  those  in  the  standard  atmosphere.  This  is 
because  the  effect  of  the  strong  temperature  inversion  at  the  tropopause 
cannot  be  properly  accounted  for  in  the  polynomial  fitting. 

Based  on  the  two  cases  presented  above,  errors  in  the  generalized 
weighting  functions  for  channels  4-7  have  negligible  effects  on  the 
retrieved  temperature  profile.  However,  a  better  polynomial  representa¬ 
tion  is  required  to  improve  the  retrieval  of  temperatures  at  the  top 
three  levels . 

3 . 7  Effects  of  Random  Errors 

Finally,  we  perform  retrieval  analyses  by  adding  random  errors  to  the 
radiance  values.  Maximum  random  errors  of  2  and  5%  were  used.  The 
resulting  temperature  errors  are  shown  in  Fig.  7.  Except  for  channel 
1,  the  addition  of  random  errors  does  not  produce  significant  errors  in 
temperature  retrievals.  This  is  particularly  evident  for  channels  4-7, 
where  peaks  of  the  weighting  function  are  in  the  lower  atmosphere.  The 
inability  to  perform  retrievals  for  channel  1  when  random  errors  are 
added  is  due  to  the  nature  of  the  broad  weighting  function.  In  this 
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case,  che  parameter  m  -  2 . 8 .  and  Ak  have  large  negative  and  positive 
values,  as  show  in  Fig.  2.  As  a  result,  small  errors  in  radiance 
values  are  amplified  and  propagate  onto  the  retrieved  Planck  intensity 


leading  to  unstable  results  for  this  channel . 


Section  4 


CONCLUSION 

In  this  report  we  have  presented  the  differential  inversion  method 
(DIM)  for  temperature  retrievals  in  terms  of  the  theory  of  radiative 
transfer.  It  is  shown  that,  in  the  Laplace  transform  space,  the  Planck 
intensity  is  directly  related  to  the  upwelling  radiance  weighted  by  the 
weighting  function.  After  expanding  the  transformed  weighting  function 
in  a  MacLaurin  series  and  performing  the  inverse  transform,  the  Planck 
intensity  in  real  space  can  be  exactly  expressed  by  a  linear  combination 
of  the  derivatives  of  radiances. 

Using  seven  HIRS  channels  in  the  15  /im  CO2  band,  we  have  performed 
numerical  analyses  for  temperature  retrievals,  including  the  curve¬ 
fitting  of  seven  radiance  values.  We  demonstrate  that  a  fifth-order 
polynomial  fitting  to  radiance  values  is  adequate  to  yield  correct 
retrieval  results  for  temperature  and  that  the  DIM  converges  to  the 
true  temperature  solution.  In  retrieval  exercises,  two  distinct 
temperature  profiles  (U.S.  Standard  and  Tropical)  were  used.  The 
retrieved  temperatures  in  the  troposphere,  corresponding  to  channels  4- 
7  are  accurate  within  1-2  K.  Addition  of  random  errors  with  a  maximum 
value  of  5%  to  the  radiance  values  does  not  introduce  instability  in 
the  inversion  exercises.  One  exception  is  channel  1,  corresponding  to 
the  center  of  the  15  /.  CO2  band,  which  has  a  broad  weighting  function. 
This  appears  to  suggest  that  the  DIM  is  particularly  useful  and 
practical  for  sharp  weighting  functions.  Our  preliminary  conclusion  is 
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Appendix  A 

DERIVATION  OF  EQ .  (2.11) 


We  rewrite  Eq.  (2.10)  in  the  form 


—  J”  e  S1T  ^  J°°  B ( 7r )  W(7r-7r)  d 7T  ^  d n 


(A.l) 


where 


r (s)  -  e  SW  R(n)  dw 


The  right-hand  side  of  Eq.  (A.l)  is  the  bilateral  Laplace  transform  of 
a  convoluted  integral  function.  Since 


-S7T  -S7T  S(ff— ?r) 

e  -  e  •  e 


(A. 2) 


We  may  rewrite  Eq.  (A.l)  as  follows: 


-  J"  e'5*  B(tt)  {  J"  w(7r__)  ^  j.  ^ 


(A. 3) 


The  teun  in  the  braces  is  the  bilateral  Laplace  transform  of  W  and  may 
be  expressed  by 


es(ff  n)  _  es>  w(y)  dy  -  w(-s) 


(A. 4) 


Substi - 


where  y  -  ir— tt  ,  dy  =  -dir,  y  -►— «  as  s  -*  ®,  and  y  -*  ®  as  ir  -»  — ® 
tuting  Eq.  (A. 4)  into  Eq .  (A. 3),  we  obtain 


r(s)  -  b ( s )  w (  s ) 


(A. 5) 


where 


b(s) 


CO  — 

e ' Sn  B (n)  dff 

—CO 


(A. 6) 


w(-s) 


CO  _ 

eS7f  W(ir)  dir 

^  —CO 


(A. 7) 


Note  that  the  integrations  of  Eqs.  (A.  6)  and  (A.  7)  are  over  the  ir- 
space.  This  is  because  the  integral  in  the  braces  of  Eq .  (A.l)  is  a 
function  of  ir  instead  of  ir.  Thus,  b(s)  and  w(-s)  are  transformed  from 
the  7r-space  to  the  s-space,  and  the  inverse  Laplace  transform  of  b(s) 
should  be  B(ir)  rather  than  B(ir) . 
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Appendix  B 


DERIVATION  OF  THE  LAPLACE  TRANSFORM  OF  THE 
GENERALIZED  WEIGHTING  FUNCTION 


On  substituting  Eq.  (2.26)  into  Eq.  (2.24),  we  obtain 


w(-s) 


m- 1 
m 


r'V) 


o 


--s+l-l  .  -1/m.  - 

p  exp  [-m  p  ]  dp 


(B.l) 


In  order  to  fit  the  integration  in  Eq.  (B.l)  to  a  form  of  the  Gamma 

1  /in  in 

function,  we  introduce  the  variable  x  -  m  p  '  ,  so  that  p  -  (x/m)  ,  and 
dp  -  m  m+^  x™  ^  dx.  Equation  (B.l)  then  becomes 


w(-s)  -  T  1(m) 


-ms  ms  -y 
x  me 


m-1 


(B.2) 


Since 


*CO 

-ms+m-1  -x  ,x  .  _  _. 

x  e  d  -  r[m(l-s) ]  ,  (B.3) 

.  0 

we  have  from  Eq.  (B.2), 

w ( — s )  -  r_1(m)  r[m(l-s) ]  mmS 


(B .  4) 
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Appendix  C 


DERIVATION  OF  EQ.  (3.18) 

From  Eq .  (3.15)  the  derivative  3E/3m  is  given  by 


!!  -  I  2[Wm(pi)  -  W(pj) 1  (C.I) 

dm  i-1  3m 

We  then  rewrite  Eq.  (2.26)  in  the  form 

wm<Pi)  “  Pi  X(m)  Y(m)  Z(m)  ,  (C.2) 

where 

X(m)  -  mm_1  ,  (C . 3) 

Y(ra)  -  r_1(m)  ,  (C.4) 

1/m 

Z(m)  -  exp  (-m  )  .  (C.5) 


Differentiation  of  Eq .  (C.2)  leads  to 


3Wm(Pi)  -  pi  (X'YZ  +  XY'Z  +  XYZ' )  ,  (C . 6) 

3m 

where 


-  X  (in  m  +  1  -  1/m)  , 

(C.7) 

-  Y  4>(m) 

(C.8) 

-  Z  p1/m  (1/m  in  pt  -  1) 

(C.9) 

On  substituting  Eqs .  (C.7)-(C.9)  into  Eq.  (C.6),  we  obtain 


wm ( P i )  Mn  m  +  1  '  1/m  ' 


3Wm(Pi)  . 


Appendix  D 


NUMERICAL  SOLUTION  OF  EQ.  (3.18) 


The  bisection  method  is  used  to  solve  Eq.  (3.18)  for  the  reason 
given  in  the  main  text.  Let  the  derivative  dE/dm  be  designed  as  f(m). 
We  start  from  an  initial  guess  -  0.1,  and  obtain  the  functii 


:ion 


value  f^  -  f(rn^).  Examination  of  the  function  dE/dm,  reveals  that 

dE/dm  >  0  for  small  m,  because  Wm(p<)  >  W(p^)  near  the  peak  and  the 

-  Vm  -  fOl  f o  1 

dominant  term  -p^  [1  -  (in  Pj.)/m]  >  0.  Thus,  if  fl  >  0,  ml  will 

be  too  small.  We  then  try  another  value,  such  that  , 

and  obtain  f^^.  If  f^  >  0,  the  process  continues  until,  for  a 

certain  ,  f^  <  0,  but  f^n  ^  >  0.  From  the  mean  value  theorem, 

the  solution  for  m  must  lie  between  rJn  ^  and  nJn^.  Subsequently,  we 

try  mtn+11  -  0.5  (m[n_1]  +  m[nl).  If  f[n+11  <  0,  we  try  m[n+2]  -  0.5 

(m^n+^  +  nJn  ^).  The  process  then  continues  until  f^m^  <  7 ,  where  7 

is  a  prescribed  convergence  criterion  and  M  the  total  number  of 

[Ml 

iterations.  The  value  m  is  then  the  solution  for  m. 
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Appendix  E 


EXPRESSIONS  FOR  THE  HIGH-ORDER  DERIVATIVES  OF  RADIANCES 

The  following  expressions  for  the  higher-order  derivatives  of 
radiances  are  derived  from  Eq.  (3.20).  Equations  (E.l)  and  (E.2)  are 
for  the  first-  and  second-order  derivatives,  respectively.  Equation 
( E . 3)  is  for  the  hth  order  derivative. 


2 

R'  -  R;  +  l  (w-wi)  R" 
1  i-1 


K 

Z  (ff-TTi)  (ff-JTj)  *  * 

k-k-2 

k>£ 


-  -  -  -  (k) 

•  (jr-irk)  R1 

(E.2) 


2  3 

Z  Z  <*-*i>  Ri’ 

i-1  j-2 
j>i 


+  •  • 


2  3 

+  Z  Z 

i-l  j-2 
J>i 


K- 1  K 

Z  Z  (w-Wj) 

i-k-2  k-k-1 
k>i 


-  -  .  ,  p(k) 

•  (*“**)  (ff-7Tk)  Rq 


(E.l) 


R"  -  2R"  +  2  Z  <*-*!>  Ri' 
1  i-l 


3  4 

+  2  Z  Z  (w_wi)  o-*j)  Ri 

i-l  j-2 

j>i 
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3  U 

+  2  z  Z 
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Z 
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>« 
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S<«  -  h!  (r{">  ♦  T  (J-Ji)  Rih+1>, 
i-1 


h+1  h+2 

+  I  I  (w-Wi)  (w-J.,) 

i-1  J-2 

j>l 


(h+2) 


i 

I 


h+1  h+2  K— 1  K 

+  I  I  •  *  •  I  I  (w-*i) 

i-1  j-2  i-k-h-1  k-k-h 

j>i  k >1 
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